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Analysis  of  a  Simple  Factorization  Algorithm 


by  (Donald  E./Knuth  iB  Luis  Trabb/Pardo 


Abstract 


The  probability  that  the  k-th  largest  prime  factor  of  a  number 


n  Is  at  most  n*  ,'  Is  shovm  to  approach  a  limit  ^  Vi  . 

Several  Interesting  properties  of  Fj^(x)  are  explored,  and  numerical 
tables  are  given.  These  results  are  applied  to  the  analysis  of  an 
algorithm  commonly  used  to  find  all  prime  factors  of  a  given  niunber. 
The  average  number  of  digits  In  the  k-th  largest  prime  factor  of  a 
random  m-dlglt  number  Is  shown  to  be  asymptotically  equivalent  to  the 
average  length  of  the  k-th  longest  cycle  In  a  permutation  on  m 


•  rC.  i  >1 1  1  y 


objects 


^  7  ^  . 
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Perhaps  the  simplest  way  to  discover  the  prime  factorization  of 
an  integer  n  is  to  try  dividing  it  by  2  ,  3  ,  h  ,  ^  .  and  to 

"cast  out"  each  factor  that  is  discovered;  we  stop  when  the  trial 
divisor  exceeds  the  square  root  of  the  remaining  unfactored  part . 

The  speed  of  this  method  obviously  depends  on  the  size  of  the 
prime  factors  of  n  .  For  example,  if  n  is  prime,  the  nimiber  of 

1/2 

trial  divisions  is  approximately  n  '  ;  but  if  n  is  a  power  of  2  , 

the  number  is  only  about  log  n  .  In  this  paper  we  shall  analyze  the 

algorithm  when  n  is  a  "random"  integer,  determining  the  approximate 

probability  that  the  number  of  trial  divisions  is  <  n  when  x  is  ^ 

a  given  number  between  0  and  1/2  .  One  of  the  results  we  shall 

35 

prove  is  that  the  number  of  trial  divisions  will  be  <  n  ,  about 
half  of  the  time. 

In  order  to  carry  out  the  analysis,  we  shall  study  the  distribution 
of  the  k-th  largest  prime  factor  of  a  random  integer.  This  problem 
is  of  independent  interest  in  number  theory,  and  for  k  >  1  it  does 
not  appear  to  have  been  studied  before.  (Wunderlich  and  Selfridge  [l4] 


gave  a  heuristic  argument  that  the  second-largest  prime  factor  will 
tend  to  be  roughly  (n^”’^^)  because  the  median  value  of 

the  largest  prime  factor  is  s«  n’^^  ;  besides  their  remark,  which 
stimulated  the  present  investigation,  the  authors  are  not  aware  of  any 
published  study  of  the  second-largest  prime  factor.  John  M.  Pollard 
[private  communication]  has  independently  investigated  the  distribution 
of  second-largest  prime  factors,  and  his  computed  values  agree  with 
those  presented  below.) 
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Section  1  of  this  paper  presents  the  factorization  algorithm 
in  detail  and  proves  its  correctness.  Quantitative  analysis  begins 

in  Section  2,  where  the  two  frequency  counts  involved  in  the  running 

* 

time  are  interpreted  in  tenns  of  the  size  of  the  largest  two  prime 
factors . 

The  distribution  of  k-th  largest  prime  factors  is  investigated 
heuristically  in  Section  somewhat  as  a  physicist  might  do  the 
analysis.  A  rigorous  derivation  of  this  distribution,  someirtiat  as 
a  mathematician  might  do  the  analysis,  is  presented  in  Section  4. 
Sections  5  and  6  continue  the  mathematical  play  by  deriving  interesting 
identities  and  asymptotic  formulas  satisfied  by  these  distributions. 
Section  7  comes  back  to  the  factorization  procedure  and  applies  the 
ideas  to  the  results  of  Sections  1  and  2,  somewhat  as  a  computer 
scientist  might  do  the  analysis. 

Section  8  discusses  the  partic\ilar  theoretical  model  used  in  these 
analyses,  and  explains  why  the  traditional  "mean  and  variance"  approach 
is  inappropriate  for  algorithms  such  as  this.  Numerical  tables  and 
empirical  confirmation  of  the  theory  appear  in  Section  9*  Finally, 
Section  10  discusses  a  rather  surprising  connection  between  prime  factors 
of  random  m-digit  Integers  and  the  cycle  lengths  of  random  permutations 
on  ra  objects. 

Although  we  shall  deal  with  a  very  simple  approach  to  factoring, 
the  results  and  methods  of  this  paper  apply  to  many  other  algorithms  as 
well.  The  paper  is  self-contained,  and  Includes  several  examples 
suitable  for  classroom  exjxjsition  of  asymptotic  methods. 
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1.  The  Algorithm. 


Here  is  the  standard  "divide  and  factor"  algorithm  which  we 
shall  analyze  in  detail.  A  proof  of  its  validity  follows  immediately 
from  the  following  invariant  assertions  governing  the  variables  used: 


n  >  2  ; 

(1.1) 

n  =  Pj^  . . 

.  m  ; 

(1.2) 

P-[^>  •  •  • » 

are  prime  nmbers; 

(1.5) 

m  >  d  ; 

(1.1) 

all  prime 

factors  of  m  are  >  d  . 

(1.5) 

Since  our  goal  is  to  analyze  a  simple  algorithm  rather  than  to  present 
it  in  optimized  form  ready  for  extensive  use,  we  shall  simply  consider 
the  following  informal  Algol -like  description: 


t  :=  0;  m  :=  n;  d  :=  2;  1 

while  d^  <  m  ^  D+l 

begin  increase  d  or  decrease  m: 

if  d  divides  m  then  D 

begin 

t  :=  t+1;  p^  :=  d;  m  :=  m/d  T-1 

end 

else  d  :=  d+1  D-T+1 

end; 

t  :=  t+1;  p^  :=  m;  m  :=  1;  d  :=  1;  1 


The  invariant  assertions  hold  after  each  line  of  this  program.  The 
expressions  in  the  rig^t-hand  column  specify  the  number  of  times  the 
operations  in  a  particular  line  will  be  performed,  where 

D  is  the  number  of  trial  divisions  performed,  (1-6) 

T  is  the  number  of  prime  factors  of  n  (counting 

multiplicity) .  (1*7) 


The  usual  refinements  of  this  algorithm,  which  avoid  a  lot  of  nonprime 
trial  divisors  by  making  d  run  through  only  values  of  the  form 
6k  + 1  ,  say,  when  d  >  3  >  have  the  effect  of  dividing  D  by  a 
constant;  so  our  analysis  of  this  simple  case  will  apply  also  with 
minor  variations  to  the  more  complicated  cases. 
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2. 


Preliminary  Analysis. 


Let  be  the  k-th  largest  prime  factor  of  n  ;  thus 

\  =  Rp+i  after  the  above  algorithm  terminates,  for  1  <  k  <  T  . 

If  n  has  less  than  k  prime  factors  (counting  multiplicities), 
let  =  1  .  We  also  let  n^  =  ®  for  convenience  in  what  follows. 

Jhe  while  loop  in  the  algorithm  can  terminate  in  three  different 
ways,  depending  on  how  we  last  encounter  it: 


Case  1,  n  <  4  .  Then  D  =  0  . 

Case  2,  n  >  4  and  the  D-th  trial  division  succeeds.  Then  the 

2 

final  trial  division  was  by  d  =  ng  ,  where  d  >  n^^  .  Since  d  is 
initially  2  and  the  statement  d  :=  d+1  is  performed  D-T+1  times, 
we  have 

D  =  ng  +  T-^  ,  ttg  >  n^^  .  (2.1) 


Case  3,  n  >  4  and  the  D-th  trial  division  fails.  Then  the 

2 

final  trial  division  was  by  d  ,  where  ng  <  d  and  d  <  n^  and 

2 

(d+l)  >  n^  .  (Note  that  if  we  set  p^  :=  1  we  have  d  >  p^  ^ 

throughout  the  while  loop.)  Thus  we  have 

D  =  fv/n^l+T-J  ,  Ug  <  n^  .  (2.2) 


In  all  three  cases  we  have  the  formula 

D  =  max(ng  ,  r/n^l ) +T  -  5  •  (2.5) 

Clearly  D  is  the  domineuit  factor  in  the  running  time,  so  most 
of  our  analysis  will  be  devoted  to  it.  However,  it  turns  out  that  the 
analysis  of  T  is  also  very  interesting;  for  large  random  n  ,  the 
number  T  of  prime  factors  can  be  regarded  as  a  normally-distributed 
random  variable  with  mean  In  In  n+1.05  and  standard  deviation 
Inn 


(see  Appendix  A) 


In  order  to  analyze  D  ,  we  shall  first  analyze  the  distributions 


of  n^  and  n^  (and  n^^  in  general) .  Let  Pj^(x,N)  be  the  number 
of  integers  n  in  the  reinge  1  <  n  <  N  such  that 


,  (5.1) 

where  x  is  any  number  >0  .  Thus  Pj^(x,N)/N  is  the  probability 
that  a  random  integer  between  1  and  N  will  have  k-th  largest 
prime  factor  <11^.  We  will  prove  that  this  probability  tends  to  a 
limiting  distribution 


Pi,(x,N) 

Ita  -  F^(x) 

N 


(5.2) 


where  Fj^(x)  has  interesting  properties  discussed  below. 

Before  we  establish  (5.2)  rigorously,  it  will  be  helpful  to  give 
a  heuristic  derivation  analogous  to  that  given  by  Karl  Dickman  [  5  ]> 
who  was  the  first  to  study  this  question  in  the  case  k  =  1  .  Let  us 


consider  Pj^(t+dt ,  N)  -  Pj^(t,N)  ,  the  number  of  n  <  N  such  that  n^^ 
lies  between  and  ,  when  dt  is  very  small.  To  count  the 

number  of  such  n  ,  we  take  all  primes  p  lying  between  and 
jjt+dt  ^  multiply  by  all  numbers  m  <  such  that  <  p  and 

mj^_^  >  p  .  Now  if  n  =  mp  we  have  n  <  and  n^^  =  p  ;  conversely 

every  n  <  N  with  nj^  between  and  will  have  the  form 

n  =  mp  where  p  and  m  have  the  stated  form.  Note  that  the  number 
of  m  <  N^”^  such  that  nij^  <  p  is  approximately  Pj^(t/(l-t)  ,  N^  ^)  , 
and  the  unwanted  subset  consisting  of  those  m  with  ^  ^  ^ 
approximately  Pj^  ^(t/(l-t)  ,  N^”^)  members.  Hence  the  number  of  m 
with  mp  <  N  and  <  P  >  P 


Pj^(t/(l-t)  ,  -  Pj^_j^(t/(l-t)  >  N^“^)  ,  ignoring  second-order  terms, 

and  we  have 

Pj^(t+dt  ,  N)  -  Pj^(t,N)  «  -  «(N^))(Pj^(t/(l-t)  ,  IT*-"*)  -  Pj^_^(t/(l-t)  ,  N^'*) 


Here  the  n  function  is  defined  as  usual, 


it(x)  =  the  number  of  primes  not  exceeding  x 


(3.3) 


(3.1^) 


According  to  the  prime  nmber  theorem  we  have  jt(x)  w  x  /  In  x  ,  hence 


.  „(N^)  «  N*dt/t  . 


(3.5) 


Plugging  this  into  the  above  formula  and  dividing  by  N  yields 


Pj^(t+dt,N)  -  Pj^(t,N)  ^^[Pj^(t/(l-t),N  "^) 

N  T  3^t 


Pl,.l(t/(l-t),N^"‘‘^) 


(3.6) 


when  N  -»  00  we  have  the  differential  equation 

Fi(t)dt  .  .  (5.7) 

Since  Pjj(O)  =  0  ,  we  may  integrate  (5*7)  to  deduce  the  formula 

*  (  (''k(3lt)-^x-l(ct))  f  •  (5-®) 


According  to  our  convention  n^  =  oo  ,  we  define 


Fq(x)  =  0  for  all  x  . 


(3.9) 


We  also  must  have 


F.  (x)  =1  forx>l,k>l  . 


(3.10) 
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Now  it  is  easy  to  see  that  (5«8),  (5.9)#  (3.10)  define  Fj^(x) 
uniquely  for  0  <  x  <  1  ,  since  we  have 

and  this  relation  defines  Fj^(x)  in  terms  of  its  values  at  points 
>  x  . 


(5.11) 
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k.  Proof  without  handwavlng. 


Our  discussion  in  the  previous  section  has  been  only  quasi-rigorous, 
but  it  shows  that  i^  the  limiting  relationship  (5-2)  holds  then 
Fj^(x)  had  better  be  the  function  defined  by  (3.8),  (5.9)>  and  (3.10). 

Now  that  we  have  a  formula  for  Fj^  ,  let  us  try  to  prove  the  limiting 
formula  (3*2). 

It  is  more  convenient  to  work  with  the  functions  defined  by 

Pj,(Q!)  =  Fj^(l/a)  ;  (k.l) 

the  above  equations  transform  into  the  somewhat  simpler  recurrence 
formulas 

Pij(Q:)  "  (Pl^(t-l) -Pj^.l(t-l))  ^  ,  for  a>l,  k>l  ;  (4.2) 

Pj^(a)  =  1  for  0<a<l,  k>l  ;  (4.3) 

Pj^(a)  =  0  for  a<0ork  =  0.  (^*4) 

Furthermore  we  let  Sj^(x,y)  be  the  set  of  positive  integers  n  <  x 
such  that  <  y  ^  and  let  Yj^(x,y)  =  llSj^(x,y)|i  be  its  csurdinality, 
so  that 

=  Yj^(N,N^)  .  (4.5) 


We  will  show  that 

Yj^(N^,N)  =  Pj^(a)N“  +  0(N“/log  n“)  , 

and  it  follows  that  a  stronger  form  of  (3.2)  is  true: 


Pk(x>N) 

N 


(4.6) 


(4.7) 
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Indeed,  we  will  prove  a  result  even  stronger  than  namely 


\(x“,x)  =  pj^(a)x“+  aj^(a)x“ /in  x“  +  O(x°‘/(log  x)^)  (4.8) 

as  X  -» 00  ,  for  all  fixed  a  >  1  ,  where  will  be  defined 

appropriately  below.  In  principle,  the  approach  we  shall  use  could 

OL 

be  extended  to  obtain  an  asymptotic  formula  for  Yj^(x  ,x)  which  is 

oj  r 

good  to  0(x  /(log  x)  )  for  any  fixed  r  ;  the  method  is  based  on 
ideas  of  N.  G.  de  Bruijn  [  1],  who  went  on  to  find  extremely  precise 
asymptotic  expansions  of  ^^(N  ,N)  in  an  elegant  way  using  Stieltjes 
integration  by  parts.  (Note;  When  k  =  1  ,  the  limiting  formula 
(3.2)  was  first  established  by  V.  Ramaswami  [11];  K.  K.  Norton  [  9] 
has  given  a  comprehensive  survey  of  the  literature  relating  to  this 
Important  special  case.) 

We  shall  use  a  strong  form  of  the  prime  niunber  theorem  due  to 
de  la  Vallee  Poussin  [  2  ] ; 

«(x)  .  L(x).0(x  ,  (4.9) 

where  C  is  a  jxjsitive  constant  and 

IW  -  4  IHT  • 

Now  to  the  proof,  which  will  be  "elementary"  except  for  our  use 
of  (4.9).  Letting  p  range  over  primes  and  n  over  positive  integers, 
we  have 
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x)  =  £  li(n  <  x“  I  =  p}ll 


X  <p  <x 


L  llCm  <  x“/p  I  <  p  and  >  p] 


X  <p  <x 


s  (\(x“/p ,  p)  -  \_i(x“/p ,  p-e)) 


a 

X  <p  <x 


where  e  is  a  small  positive  number  and  YQ(x,y)  =  0  .  The  key  idea 
in  o\ir  derivation  will  be  to  replace  £  /p  >  p)  by 

X 


X  <p  <x 


a 


a 


J  '4'j^(x^/y ,  y)dy/(ln  y)  ,  using  the  ''density'*  function  for  primes 


suggested  by  (4.10).  To  justify  this,  we  have 

^x<p<x  J  ^ 


a 


x<p<x“  n€S.(x®/p,p) 


£  X  - ; 


E 


Qt-l 


1  <n  <x 

OCi 

nj^<x  /n 


r 


V 


A 


\<P<x“/n 

x<p  y 


*  '  n€S,(x“/y,y) 


X  /n  J 

r  '  J2_ 

J  In  y 

max(nj^,x) 


1 1  ^ 

In  y 


J 


£  ^x(x°^/n)  -  «(max(nj^,x))  +  0(1)  -  L(x^/n)  +  L(max(nj^,x)  )^ 


1  <n  <x 
i\<x  /n 
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\<X  /n 

=  O(x^(log  ^ )  .  (^*11) 


A  similar  estimate  applies  to  Z/  a  \  1^^  have 


X  <p  <x 


a 


\(x  ,x)  =  X 


»/('.(7-)-'«(f .)) 


JSL- 

In  y 


-r 


OCX. 


a 


\^(lOg  x) 


0 


(1|.12) 


as  X  -•  ®  ,  for  all  fixed  r  >  0  .  This  is  the  formula  we  shall  use 
for  Q!  >  1  ;  for  0  <  a  <  1  we  have  'l'j^(x'^,x)  =  •  (The 

brackets  L  J  in  the  latter  formula  turn  out  to  be  important,  since 

the  integral  (1+.12)  is  sensitive  to  0(1)  terms  in  the  vicinity  of 

a  . 
y  =  X  .) 

Our  proof  of  (^+.8)  is  by  induction  on  k  ,  and  for  fixed  k  by 
induction  on  reel  .  Actually  the  first  case  k  =  l,  fo!!  =  2  seems 
to  be  the  hardest ;  when  1  <  a  <  2  we  have 


15 


where  (x)  denotes  x-Lxj  .  The  remaining  integral  is 


=  lim  ((In  n)  -  (H^-1))  =  I-7  , 

n  —<0 


(4.14) 


where  7  is  Euler's  constant. 

Now  suppose  we  have  proved  that 


Y^(x“  X)  .  x“  P^(a)  +  (1-r)  ^ 


(‘‘•15) 


for  1  <  a  <  m  ,  where  the  bounding  constants  implied  by  the  0  '  s  depend 
on  m  but  not  on  x  or  a  .  The  discussion  in  the  previous  paragraph 
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establishes  (^^-15)  for  ra  =  2  .  We  can  extend  it  to  the  next  case  by 
analyzing  its  value  for  m  <  a  <  nH-l  ; 


““  4  ‘  1  x)  )  f 

4  "( 


a/t 
fa/tjiog  X 


.(log  x)‘ 


(lv.l6) 


by  substituting  x^/^  for  y  and  inserting  (I.15).  Continuing,  we 
get 


.  x^-x®  f  Pj(t-l)  f  *  ^a/t  M 


^  a 

X 


V(log  x>‘ 


a/2 


a 


f  ■>.  /I  \  a  a-l 

a  ,„A  .  a  luldu  (l-7>x  c  dt 

=  X  P  (a)  +  X  r  g  ;  a,  '  -  “S’  Ji  ”1'  ’  “ 

1  u  ln(x  /u)  In  X  1 


\^(log  x)‘ 


0(7^1  , 

In  X  \.(l08  y 


(1^.17) 


since 


a/2 

f*  (u}du 

1  ln(x*^/u) 


In  X 


a 


as  in  (if-. 13)  and  (4.lU),  and 


0 


x°^  du 

V  j 


(U.18) 


(4.19) 


with  bounding  consteints  depending  only  on  m  .  This  establishes  (4.15) 
for  all  m  ,  by  induction. 

We  have  proved  (4.8)  for  k  =  1  ,  with 


o^(a)  =  (l-7)p^(a-l) 

For  larger  k  ,  a  similar  but  simpler  derivation  applies:  Assuming  that 


Vx“x)  =  P^(a)  .  ^ 


for  1  <  a  <  ra  (cf.  (4.15)),  we  extend  this  to  m  <  a  <  nH-1  by 


(4.20) 


\(x“x) 


(X 


+  0 


Jlog  x)' 


a 

=  X 


a 


dt 


j“  (p  (t-l)  -  P^.i(t-l))  ^  I  Y  0  (_J^ 
In  2  ^  }  Vice  x)‘ 


;  (1^.21) 


i6 


the  desired  relation  follows  for  k  >  2  provided  that  we  define 


=  0  for  a  <  2  . 

It  follows  that 

=  (l-7)(pj^(Q!-l)  -  Pj^_j^(Q!-l)) 
for  all  k  >  1  • 


(4.22) 

(4.23) 

(4.24) 
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5-  Identities  satisfied  by  pj^  . 


The  functions  defined  by  {h.2),  (h.k)  possess  many 

rather  surprising  properties,  and  we  shall  examine  some  of  them  in 
this  section. 

Our  first  goal  is  to  express  the  pj^  in  terms  of  the  poly  logarithm 


functions 

Lj^  ,  defined  by 

Lo(a) 

=  0  for  a  <  0  , 

L^Ca)  =1  for  a  >  0  ; 

(5.1) 

(X 

-  4  f  • 

(5.2) 

a 

Thus  L,  (a)  =  In  a  for  a  >  1  ,  and  Lp(a)  =  f  ln(t-l)dt/t  for 
a  >2  ,  etc.;  it  is  not  difficult  to  verify  that  is  l/kl 

where  1  <  x^, . . .  ,x^  <  a  and  1  >  ^  i  /  d  •  In 

particular,  =  0  for  a  <  k  . 

By  iterating  the  recvirrence  for  pj^  we  find 

l-Pj^(a)  =  -L^{a)  +  L^(a)  -  Lj^(a)  +  L^(a)  -  ...  ,  (5.5) 

1  -  p2(a)  =  LgCa)  -  2Lj(a)  +  5L^(a)  -  4L^(a)  +  . . .  ,  (5.4) 

for  a  >  0  ,  and  in  general 
=  n>0  ^ 

These  infinite  suras  are  actually  finite  for  any  peurticular  value  of  a  . 

Now  let  us  examine  several  auxiliary  functions: 


.Xj^)  over  all  points  x^,  ...,Xj^ 


times  the  integral  of  (dae2^. .  .dXj^)/( 
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for  p  >  a  or  p  <  0  ; 


(5.6) 


Si.(a,p)  =J 


a  pj^(t-l)dt 


p-t 


Sj,(Q!)  =  Sj^(a,QH-l)  ; 


(5.7) 


a  p  (t-1) 


=  f 


dt  ; 


(5.8) 


(5.9) 


-(x)  =  I  Pk(t)e“'^^dt 


>  X  ^  0  * 


(5.1a) 


(This  is  a  different  function  fro®  that  in  Section  k.)  It  follows 

immediately  from  the  definition  P}j(of)  =  1  -  ‘^j^(^)  +  l(<^)  that 

<^1^(0!)  =  k  -  p^(a)  -  . . .  -  pj^(o;)  .  (5. 11) 

Integration  by  parts  enables  us  to  evaluate  Ij^(a)  as  follows; 

a  «  P],(t)dt 

ik(«)  -  w«)  =  -  Pk(^)  ^  J 


=  -  Pij(a)  In(OH-l)  +  aj^(offl)  . 


(5.12) 


Thus  in  particular  we  have 


I^(a)  =  -  p^(a)  In(cH-l)  +  1  -  pj^(oH-l)  , 


(5.15) 


Ipicc)  =  -  Pi(Q!)  lii(oH-l)  -  p2(a)  ln(QH-l)  +  5“2pj_(QH-l)  -  pgCoH-l)  ,  (5.1k) 


etc.  A  somewhat  surprising  consequence  of  this  relation  is  that 
Ij^(oo)  =  k(k4-l)/2  ,  while  orj^(®)  =  k  ;  in  particular,  Ij^(®)  =  o^(<x>) 
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yields 


Integration  by  parts  applied  to  Sj^(a,p) 


P-t 


a 


a 


ap  (a)  QH-l  p  (t-l)dt 

P  J  - - 2"  • 

•'l  (p+l-t)^ 


(5.15) 


Differentiating  the  integral  which  defines  Sj^(a)  =  Sj^(a,(>H)  with 
respect  to  a  leads  to  a  formula  which  can  be  combined  with  this  one; 


a  p  (t-l)dt 


=  Pij(a-l)  -  i  ((Q!-l)p^(a-l)  +  Sj^(a-l) -Sj^_^(o£-l)) 


=  5  (pj^(Q!-l)  +  Sj^.3^(a.l)  -Sj^(a-l))  .  (5.16) 

Now  we  are  ready  to  prove  an  Important  relation  which  expresses  p^^^ 
in  terras  of  pj^  and  pj^  ^  ; 


Lenma. 

i  '  for  k>l  .  (5.17) 

Proof.  Since  ~  Pk^°^^  ~  ^  ® 

0  <  a  <  1  ,  the  result  holds  for  Foci  =  1  ;  we  will  show  that  the 
derivatives  agree,  by  induction  on  fal  •  Since 
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(QH-I)p^^(an)  =  Pk(«)  -  Pfcfi(«)  =  (Sj^_3^(a) -Sj^(a))/k 

(Qn-1)  p^(an)  =  Pj^_i(Q!)  -  pj^(a)  , 

(QH-l)s^(aH)  =  pj^(a)  +  Sj^_^(Q!)  -  Sj^(a)  , 

(Qffl)S^_^(QH-l)  =  Pk-l(“)  +  \.2(“)  » 

the  desired  result  is  equivalent  to 

¥  fk(“)  =  ¥  *  i  ("k-if^)  • 

For  k  =  1  this  is  obvious,  otherwise  it  holds  by  Induction.  □ 

By  iterating  the  recurrence  in  the  lemma,  it  follows  that 

Finally  let  us  consider  the  functions  ej^(x)  defined  in  (5. 10). 
Somewhat  surprisingly,  these  can  actually  be  expressed  in  closed  form: 

Theorem.  ej^(x)  =  —  ^1  +  +  ...  +  ^  ,  where 

E(x)  =  Ej^(x)  is  the  exponential  integral  function 
00  00 

E(x)  =  J  e'*  dt/t  =  f  e"^  dt/t 
X  1 


(5. 18) 


(5.19) 
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Proof.  Once  again  we  integrate  by  parts: 


®  p.  (t-1)  -  (t-1)  ./t.i'i, 


dt 


-  e*  /  t  e""*^*  d  p  (t) 
0 


/  Pj^(t)(e"  *-txe"'^^)dt 
0  ^ 


=  e  (ej^(x)  +  X  e^(x))  . 

E(x) 


If  we  let  fjj(x)  =  xe  ^  have  therefore 

^L(^)  =  (x)  +xe‘(x)  -e‘*e  (x)) 


-X 


and  it  follows  by  induction  on  k  that 
f  (x)  -  C  + 

In  order  to  evaluate  C  ,  we  integrate  by  parts  in  the  opposite  direction: 


X  ej^(x)  = 


-  /  pJt)d(e-**)  = 

r\  ^ 


;V**4p  (t) 
0  0  ^ 


'  1  -  2 


dt 


=  1 


Hence  C  =  lira  x  e,  (x)  =  lira  f,  (x)  =  1  .  □ 

x-*oo  k^  '  X k'  ' 
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6.  Asymptotic  formulas. 

In  this  section  we  shall  study  the  asymptotic  behavior  of  Pj^(Q!) 
for  large  a  .  Our  starting  point  is  a  simple  proof  that  pj^(Q!)  is 
exponentially  small:  Let  us  write  p(a)  for  p^(o!)  .  Then  since 

a  a 

1+J  p(t-l)dt  =  J  p(t-l)dt 


a  a 

=  -  tp(t)  I  +  r  p(t)dt 

1  1 

QH-l 

=  l-ap(oi)  +  J  p(t-l)dt 


(6.1) 


we  have 


Qifl 

J  p(t-l)dt  =  ap(a)  .  (6.2) 

a 

It  follows  immediately  that  Q!p(ct)  <  p(q:-1)  for  all  a  >  1  ,  hence 
by  induction 

p(n)  <  l/nl  (6.3) 

for  all  integers  n  >  1  .  Considerably  more  precise  fomulas  have  been 
obtained  by  de  Bruijn  [  1  ]  and  others,  and  niunerical  results  have  been 
tabulated  by  Mitchell  [8]  and  by  van  de  Lune  and  Wattel  1 15];  but 
(6.3)  suffices  for  our  purposes  in  this  section. 

The  rapid  decrease  of  p^(ct)  simplifies  the  numerical  evaluation 
of  integrals  and  it  also  leads  to  a  simple  treatment  of  the  asymptotic 
behavior  of  P2(Q!)  : 
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Theorem .  For  all  fixed  r  >  1  we  have 

p  (a)  =  A^— +  V 

\  a  a  a  J 


as  a  -•  00  ,  where 


(6.4) 


A  =  e'  w  1.78107  24179  90197  98524  , 


(6.5) 


and  the  coefficients  Cj^  are  defined  by 

z 


£  z^Cj^kl  =  exp^J*  (e^-l)  dt/t  ^  =  exp^  S  z^/k-k'.^  . 


k>0 


(6.6) 


Thus  <Cq,c^,C2,...>  =  ^1,  1,|  ,  ^  ^  ,  y  ,  ^  , 

proving  the  theorem,  we  note  that 
(6.6)  implies  the  recurrence  formula 


n 


1 

n 


E 

l<k<n 
n-1 


("Jv 


n  >  1 


(6.7) 


Therefore  c  >  c  „  for  n  >  2  ,  and  c,_  ^  >  nl  ;  the  infinite 
n  2  n-2  —  '  2n^•l  ’ 

series  E  diverges  for  all  a  .  In  other  words,  (6.4)  is 

strictly  an  asymptotic  formula. 


Proof. 


P2(«)  = 


From  the  lemma  in  the  previous  section  we  have 

=  p(o).j“ 

-  p(c.)  ■>  I  pCt-ijdtf  j: .  ■>  ■  ■  ■  -f  ' 

1  \cc  a  cr^-^  ar'^(Qffl-t) 


a-1 


J  p(t)t^  dt/ot^^  +  0(a"^"^) 

0 


E 

0  <k<r 


(6.8) 


a 


since  f  p(t-l)  (t-1)^^  dt/(Qrt-l-t)  <  f  p(t-l)  (t-1)^^  dt  <  ® 
1  1 


Furthermore  we  have 

m 

k 


-ia 

/  P(t)t^dt  =  o(  /  e‘'‘^t^dt  1  =  ol  e  ^ 


a-1 


a-1 


(6.9) 


as  a  -*  ®  ,  by  making  very  crude  estimates  not  even  as  powerful  as 
(6.5),  so  we  can  integrate  to  ®  in  (6.8): 

p_(a)  =  -0  +  -i  +  . . .  +  ^  +  o(a’^“^)  , 


(6.10) 


a  a 


a 


where 


\  =  /  p(t)  t^dt  . 

0 

It  remains  to  evaluate  the  .  We  have 


(6.11) 


S  =  /pWe-'^dt  .  e,(x)  =  (6.12) 

k>0  ^  0  ^ 


by  the  theorem  of  Section  5;  and  it  is  well  known  that 


-  E(x)  -  In  X  =  7  +  £  (-'x)^/k.k'.  .  (6.15) 

k>l 

(See,  for  example,  [7,  exercise  5.2.2-U5].)  This  combines  with  (6.12) 
and  (6.6)  to  prove  that  a^  =  e  Cj^  .  □ 
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The  coefficients  Cj^  have  the  curious  property  that 


(6.14) 


is  also  an  asymptotic  expansion  of  ,  but  not  as  accurate  when  truncated. 
Another  series, 


P2(«) 


=l-°0  .  . 

- 2  *  - T 

2(a-l)  3{ci-iy 


is,  in  turn,  more  accurate  than  (6.4) .  These  series  are  obtainable  from  one 
another  using  the  relation  P2(ci!)  =  -(Q6fl)p^(ot<-l)  +  p^(a)  . 

For  k  >  5  i  we  shall  content  ourselves  with  establishing  the  leading 
term  in  the  asymptotic  exx>ansion  of  pj^  ,  namely 


(6.15) 


[Appendix  B  contains  an  asymptotic  expansion  of  pj  .  ]  Consider  first 


«  ■  1  (e  •  «(j)] 


dt 


Offl-t 


and  note  that 


a 


dt 


a 


4  t^(cw-l-t) 


cw-l 


dt  p 

dt 

\  2  In  a 

h 

OH-l-t 

1  ■  QH-l 

dt  1 

a 

dt 

4 

t(ofH-t) 

.  i  1  +  2 

In  a 

(6.17) 


Hence  82(0!)  =  2  A a"^  In  a  +  0 (a"^)  ,  and  pj(Q:)  =  Aa^^ln  a+0(a*^) 
by  (5.18).  In  order  to  use  this  approach  for  larger  k  ,  we  note  that, 
when  k  >  1  , 


f’'^  (In  t)^  dt 
I  t(Qf4-i-t) 


1  [»  (In  t) 

Offl  i  t 


^  (In  t)^  dt  1  ^  (In  t)^  dt 

t  OW-1  QH-l-t 


I  (111  ar- 

■(ofiy  (k+i)  Qffi 


k  (In  ln(cefl-t)dt 


1  (In  ar 
(ofl)  (fcfl) 


(In  ln(l  -  ^)dt 


Now  In(l-x)  =  -  X  f (x)  ,  where  f  is  a  function  satisfying 

1  <  f(afi)  < 


(6.18) 


when  1  <  t  <  a  ,  hence 


4  ln(  1  -  ^  )dt  =  ^  j“  (On  0 


(In  a)dt 


=  0(ln  a)^  . 


(6.19) 


We  have  proved  ■chat 


(In  t)  dt 
t(CH'l-t) 


(6.20) 


for  all  k  >0  .  Using  (5.18),  formula  (6.15)  now  follows  by  induction, 
together  with 


\(C,)  . 


Ak(ln 


1):  a  V  ®  / 


(6.21) 
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7-  Application  to  factoring. 

The  distributions  Fj^(x)  =  pj^(l/x)  can  be  used  to  estimate  the 

running  time  of  various  algorithms  for  factorization.  For  example, 

Pollard's  important  new  Monte  Carlo  method  [10]  takes  about  steps, 

where  is  the  second-largest  prime  factor  of  n  ,  so  we  can  use  a 

table  of  Fg  to  state  that  Pollard's  method  will  complete  the 

factorization  in  0(n‘^*^^)  steps  at  most,  about  half  of  the  time. 

For  the  simple  algorithm  of  Section  1,  we  need  to  analyze  the 

distribution  of  max(n2,“^^)  ,  and  this  does  not  appear  to  be 

expressible  directly  as  an  algebraic  function  of  the  Fj^  .  However, 

we  can  readily  carry  out  the  analysis  by  using  the  techniques  above. 

Let  G(x)  be  the  limiting  probability  that  max(n2  ,  Vn^ )  <  N*  , 

when  n  is  a  random  integer  between  1  and  N  .  Then 

G(x)  =  Fj^(x)  +  G^(x)  =  FgCx)  -G2(x)  ,  where  G^(x)  is  the  probability 

that  and  Ug  <  ,  emd  Gg(x)  is  the  probability  that 

2x  jc 

n^  >  N  and  Hg  <  w  .  Arguing  as  above,  we  find 
2x  r  ^  \  2x  /'T  4.  \ 

(7.1) 


(7.2) 


2x 


t  i-t  ) 

I  t  pU  J 

2/a 

r  p((l-t)a) 
ya 

dt 

t 

p(u-l)du 
-  4.3^  c»l-u 

t  ^li^l-t  ) 

4  t  pU  j  ^ 

(7.5) 


(7.1^) 


(note  that  Cii(i)t(52(i)  =  S^(a)  =  F2(i)-F^(i)  ,  in 
agreement  with  the  lemma  of  Section  5*)  It  is  clear  from  our  asymptotic 
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results  that  G^(l/a)  decreases  exponentially  for  large  a  ,  hence 
it  is  numerically  better  to  use  the  formula  G(x)  =  Fj^(x)  +  Gj^(x) 
than  to  use  F2(x)  -G2(x)  ;  furthermore  the  integration  is  over  a 
limited  range.  On  the  other  hand  for  2  <  a  <  5  it  is  most 

ln(a/2)  in  this  range. 


convenient  to  use 


G2  since  Gg, 


(i)  = 
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8.  Remarks  about  the  model. 


The  probability  considerations  above  are  for  random  n  between 
1  and  N  ,  and  for  relations  such  as  <  N*  ;  but  from  an  intuitive 
standpoint  we  might  rather  ask  for  the  probability  of  a  relation  such 

X 

as  <  n  ,  without  considering  N  .  Actually  it  is  easy  to  convert 
from  one  model  to  the  other,  since  most  numbers  between  1  and  N  are 
large. 

More  precisely,  consider  how  many  numbers  n  between  ^  N  and  N 
have  n^  <  ;  this  is  Pj^(x,N)  -  Pj^^  ^  >  I  )  =1  N*Fj^(x)  +  0(N  /  log  N)  j 

since  Pj^(x,N)  =  N*Fj^(x)  +  0(N  /  log  N)  .  Furthermore,  consider  how 
many  of  these  n  have  n^  <  n^  <  :  The  latter  relation  implies 

N^-nic>(|N)  2/log  N  ^  ^ 

Fj^(x)  +  0(1  /  log  N)  ,  since  Fj^  is  differentiable;  so  the  number  of  such 
n  is  at  most  Pj^(x,N)  -  Pj^(x  -  log  2/log  N  ,  N)  =  0(N  /  log  N)  .  (The 
constant  implied  by  the  0  in  (U.7)  will  be  independent  of  x  in  a 
bounded  region  about  x  .) 

We  have  shown  that  Fj^(x)  +  0(l/  log  N)  of  all  n  between  ^  N 
and  N  satisfy  n^^  <  n  .  Therefore  if  Qj^(x,N)  denotes  the  total 

X 

number  of  n  <  N  such  that  n^^  <  n  ,  we  have 

Q  (x,N)  =  E  -^n(f(x)  +  o/^ - V 

^  l<j<log2logN  2^  \log(N/2^)yy  V  log  N  y 

by  dividing  the  range  N/log  N  <  n  <  N  into  logg  log  N  parts . 


(8.1) 
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It  is  customary  to  define  the  "probability"  of  a  statement  S(n) 
about  the  positive  integer  n  by  the  formula 

Pr(S(n))  =  lim  ^  (nmber  of  n  <  N  such  that  S(n)  is  true)  ,  (8*2) 

N-oo  ^ 

when  this  limit  exists.  Thus,  we  can  state  well-known  facts  such  as  the 
following:  Pr(n  is  even)  =  ^  ;  Pr(n  is  prime)  =  0  ; 

Pr(n  is  squarefree)  =  6/ jt  .  Equation  (8.1)  now  yields  another  result 
of  this  type: 

^ 

for  all  fixed  x  . 

Another  important  observation  should  also  be  made  about  the  theoretical 

model  we  have  used  to  study  the  factorization  algorithm  in  this  paper: 

We  have  stated  our  results  in  terms  of  the  probability  that  the  running 

time  is  <  (or,  if  we  prefer,  n*  ) ;  this  contrasts  with  the  customary 

approach  to  the  study  of  average  running  time,  which  derives  mean  values 

and  the  standard  deviation.  The  reason  for  abandoning  the  traditional 

approach  is  that  the  mean  and  standeurd  deviation  are  particularly 

uninformative  for  this  algorithm.  This  phenomenon  is  apparent  when  we 

consider  that  the  mean  running  tirae^over  all  n  <  N  will  be  relatively 

near  the  worst  case  n  ,  but  in  more  than  70  per  cent  of  all  cases  the 

0.4 

actual  running  time  will  be  less  than  n  *  . 

In  order  to  understand  this  rather  anomalous  situation  more  fully,  let 
us  calculate  the  asymptotic  mean  and  standard  deviation  of  the  largest 
prime  factor  n^^  ,  when  all  integers  1  <  n  <  N  are  considered  equally 
likely.  Let  $(t)  be  the  probability  that  <  t  ,  when  n  is  in  this 
range.  Then  the  derivation  of  Eq.  (4.13)  allows  us  to  conclude  that 
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$(t) 


(8.4) 


for  <  t  <  N  . 

We  shall  now  calculate  the  asymptotic  behavior  of  the  k-th  moment 

k 

of  this  distribution,  namely  the  asyn5>totic  expected  value  of  n^  . 
[Incidentally,  our  derivation  provides  a  good  example  of  the  use  of 
Stieltjes  integration.]  The  k-th  moment  is 
N 


E(np  =  f  t'^di(t) 


and  since  the  integral  from  1  +  /n  is  f  d$(t)^  =  0(N^/^) 

it  can  safely  be  ignored.  We  are  left  with 


(8.5) 


(8.6) 


by  replacing  t  by  N/v  in  the  second  integral.  [The  0  estimate  here 

b 

is  justified  by  the  following  general  lemma:  Let  f  f(t)  dg(t)  and 
b 

f  f(t)  dh(t)  exist,  where  h(t)  =  0(g(t))  ,  and  where  both  f  and  g 
a 

are  positive  monotone  functions  on  [a,b]  .  Then  it  is  easy  to  see  that 

b  /  b 

/  f(t)  dO(g(t))  =  0(f(a)g(a))  +  0(f(b)g(b))  +  o(  f  f(t)  dg(t) 
a  V  a 


if  we  integrate  by  parts  twice.]  The  first  integral  in  (8.6)  is 
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N 


; 


In  t 


=  n"; 


dv 


1  N  -  In  v) 


In  N 


/n 


f "  +  r  _ 

J-  k+l  «,  k+l/T  „  T  V 

1  V  1  V  (InN-lnv) 


y? 


In  V  dv 


k  In  N 


\(log  N)‘ 


) 


yN 


Jc  f*  kf2 

The  second  integral  is  -a /in  N  times  J  (v}dv/v  >  which  is  within 
0(N-'’^^>/=)  of 


J 


[vldv 

,  k+2 

1  V 


,  T  J.  k+2 

d  >1  J  V 


iwi 


irr  -  IIT  (C(w-i)-i)  =  i  -  . 


k(k+l)  "  k+l 


Thus  we  have  shown  that 


e(A  ,  gi^i)  jii  +  of 
^l’  l""  V(X=gN)V 


(8.7) 


It  follows  that  the  mean  value  of  n^^  is  asymptotically 
(it^/l2)N/ln  N  ,  and  the  standard  deviation  is  (C(3)/5)^^^N/Vln  N  ,  to 


within  a  factor  of  1+0(1/  log  N)  .  In  particular,  the  ratio 
standard  deviation 

'  -4  CO 

mean 

as  N  -♦  »  ;  this  result  danonstrates  the  unsuitability  of  a  traditional 
"mean  and  variance"  approach  to  the  analysis  of  such  algorithms. 


(8.8) 
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9-  Numerical  results. 


The  differential-difference  equations  for  are  conveniently 

suited  to  numerical  integration.  For  example,  given  internal  arrays 
containing  pj^(m  +  k/n)  ,  p2(m  +  k/n)  ,  and  p^(m  +  k/n)  for 
0  <  k  <  n<-t  ,  where  m  is  some  fixed  integer  and  B  =  l/n  is  the 
step  size  and  t  depends  on  the  method  of  integration,  one  pass  over 
these  arrays  serves  to  increase  m  by  1  .  When  m  reaches  a 
suitably  large  value,  the  asymptotic  formulas  derived  above  provide 
an  excellent  check  on  the  accuracy  of  the  calculations.  Another 
excellent  check  comes  ftom  the  formula 
00 

e’'  =  J  p(t)dt  =  p(l)  +  2p(2)  +  5p(5)+ ...  ;  (9.1) 

0 

cf.  (6.2),  (6.5),  and  (6.11).  (Incidentally,  identity  (9.I)  appears  to 

• 

be  new;  it  was  discovered  empirically,  after  noticing  that  the  results 
of  nvimerical  integration  seemed  to  resemble  a  "familiar"  constant.  This 

y 

particular  constant  came  as  a  surprise,  since  e  usually  occurs  only 
in  connection  with  infinite  products.  After  the  proof  of  (9.I)  was 
found,  the  theorem  in  Section  5  above  followed  rather  quickly.  Thus, 
numerical  results  indeed  suggest  theorems.) 

The  following  table  gives  representative  values  of  p^  ,  p^  ,  p^ 


and  G  to  12D  : 


a 

Pl(a) 

P2(«) 

PjCa) 

G(l/a) 

1.0 

1.000000  000000 

1.000000  000000 

1.000000  000000 

1.000000 

000000 

1.5 

.59U534  891892 

1.000000  000000 

1.000000  000000 

1.000000 

000000 

2.0 

.306852  81941IO 

1.000000  000000 

1.000000  000000 

1.000000 

000000 

2.5 

.130319  561832 

.953389  706294 

1.000000  000000 

.750246 

154979 

5.0 

.048608  388291 

.852779  323041 

1.000000  000000 

.447314 

214952 

3.5 

•016229  593243 

.735481  165219 

.997526  275042 

.225819 

495955 

Iv.O 

.004910  925648 

.625681  059959 

.985113  653272 

.096399 

005935 

1I.5 

.001370  1177^+1 

.533652  572034 

.960975  011157 

.036573 

065077 

5.0 

.000354  724700 

.463222  186987 

.927859  655628 

.012413 

482748 

6.0 

.000019  649696 

.365217  751694 

.851107  195638 

.001092 

266742 

7.0 

.000000  874567 

.301786  010308 

.777229  329492 

.000071 

391673 

8.0 

.000000  032321 

.257455  710831 

.712844  794121 

.000005 

662651 

9.0 

.000000  001016 

.224592  162720 

.657959  58195*^ 

.000000 

155284 

10.0 

.000000  000028 

.199248  208994 

.611115  9975k) 

.000000 

005585 

12.0 

.000000  000000 

.162658  856655 

.555865  615616 

.000000 

000004 

lU.O 

.000000  000000 

.137457  568144 

.478221  749442 

.000000 

000000 

16.0 

.000000  000000 

.119016  455055 

.452642  865552 

.000000 

000000 

18.0 

.000000  000000 

.104958  755569 

.395653  753569 

.000000 

000000 

20.0 

.000000  000000 

•095875  845625 

.364991  546696 

.000000 

000000 

25.0 

.000000  000000 

.074277  805044 

.507069  057805 

.000000 

000000 

30.0 

.000000  000000 

.061455  756517 

.266170  912880 

.000000 

000000 

ItO.O 

.000000  000000 

.045685  815582 

.211858  770558 

.000000 

000000 

50.0 

.000000  000000 

.056556  095670 

.177085  969207 

.000000 

000000 

60.0 

.000000  000000 

.030192  05?r52 

.152778  425205 

.000000 

000000 
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Figure  2.  Empirical  distribution  functions  corresponding  to  Figure  1, 
based  on  the  factors  of  the  largest  100  10-digit  numbers . 


(In  1950,  Diclonan  published  8d  values  of  integer 

a  <  8  ;  his  figures  were  correct  except  that  pj^(7)  was  given  as 

"  .0000  0088  ”.) 

Figure  1  shows  these  distributions  graphically,  and  illustrates 
the  fact  that  F^(0)  =  G' (O)  =  F^(  |  |  ^  =  0  ,  F^(0)  =  A  , 

=2  ,  F^(l)  =  1  ,  F^(0)  =  00  .  Although  the  graphs  of  F^^  , 

F^  )  and  F^  are  qualitatively  different,  the  graphs  of  Fj^  for 
k  >  4  will  resemble  that  of  F^  (but  they  will  rise  ever  more 
steeply) . 

The  following  table  shows  percentage  points  of  the  distributions 
Fg  ,  Fj  ;  for  example,  the  probability  is  only  10  percent  that 

,18616 

n 


^1^ 
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p 

Fg^Cp) 

Fj^Cp) 

.01 

.26971I 

.00558 

.00068 

.02 

.29341 

.OHIO 

.00149 

.03 

.31004 

.01656 

.00239 

.Oil 

.32341 

.02196 

.00334 

.05 

.33483 

.02730 

.00435 

.10 

.37851 

.05308 

.00995 

.15 

.41288 

.07741 

.01629 

.20 

.44304 

.10033 

.02327 

.25 

.47068 

.12191 

.03079 

.30 

. 49656 

.14216 

.03882 

,ko 

.54881 

.17892 

.05636 

.50 

.60653 

.21172 

.07584 

.60 

.67032 

.24267 

.09745 

.70 

.74082 

.27437 

.12165 

.75 

.77880 

.29153 

.13506 

.80 

.81873 

.31035 

.14972 

.85 

.86071 

.33201 

.16627 

.90 

.90484 

.35899 

.18616 

.95 

.95123 

.39672 

.21377 

.96 

.96079 

.40681 

.22141 

.97 

.97045 

.41850 

.25054 

.98 

.98020 

.43268 

.24224 

.99 

.99005 

.45169 

.25954 

1.00 

1.00000 

.50000 

.55335 

Empirical  confirmation  of  the  theory  is  illustrated  in  Figure  2, 
which  shows  exact  empirical  distribution  functions  corresponding  to 
Figure  1  for  the  100  nxanbers  n  =  10^^ -m  ,  1  <  m  <  100  .  As  expected, 
the  deviation  frcan  is  most  pronounced  for  k  =  1  and  x  >  ^  , 

but  the  deviations  are  not  severe.  This  set  of  numbers  contains  three 
primes  (10^*^  -  55  ,  10^^  -  57  >  10^  "71)  >  and  ten  products  of  two 
primes.  The  smallest  values  of  n^^  occurred  for 
10^° -100  =  157-101*75*11-5^*3^*2^  ,  10^° -  64  =  465-451-29*5^'2^  ; 

the  largest  values  of  n^  occurred  for  10^-69  =  456767 '21895  , 

10^^  -  22  =  85021*19605  •5*2  ;  the  largest  values  of  n^  occurred  for 
10^° -51  =  88501'421‘269  ,  10^° -75  =  15879- 559 -225 *5^  •  The 

smallest  values  of  max(Vn^  ,  n^)  occiirred  for 

10^°-100  =  157’101-75*11'5^-5^‘2^  ,  10^°-25  =  2857'115-59'7*5^-5 
(so  these  would  be  the  easiest  numbers  in  the  given  range  to  factor  by 
the  simple  algorithm);  the  smallest  values  of  n^^  for  which 
'/n^  >  Ug  occurred  for  10^^  -  66  =  59^17  *105  *45 '19  *2  , 

10^-68  =  77201 -55 -47 -15  *2^  . 

In  Dickman's  original  paper  he  calculated  the  "average"  value 
of  X  such  that  n^  =  n  ,  namely  the  exi>ected  value  of  log  log  n  . 

This  equals 

1  00  CO 

D,  =  r  xdF,(x)  =  -  J  p'(t)  dt/t  =  J  p(t-l)  dt/t^  (9.2) 

0  1  1 

and  by  Eq.  (5'1*<^)  we  also  have 

I 

00  00 

r  p(t-i)  dt/t^  =  -  s,  (00,-1)  =  J  p(t-i)  dt/(t+i)  .  (9.5)  I 

1-^1  * 

'  j 

j 

I 
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In  a  similar  way  we  cein  determine  the  expected  value  of  log  /  log  n  , 
a  nimher  vrtiich  can  be  expressed  in  several  ways,  namely 

1  00  ® 

\  =  [  (pj^(t-l) -Pj^_3^(t-1))  dt/t^  =  1  -  J^Pi,(t)  dt/t^ 


=  (Pl,(^-1)  -2Pk.i(t-l)+Pl,_2(^-l))  dt/Ct+l)  .  (9.4) 


Nmerical  evaluation  (using  the  asymptotic  formulas  for  p^  and  ) 
gives 


=  .62432  99885  ; 

Dg  =  .20958  08743  ; 

Dj  =  .08831  60989  . 


(9.5) 

(9.6) 
(9-7) 


(Dickraein's  value  for  was  .624529998  .  Note  that  Dg  is  not  equal 
to  D^(l-D^)  ,  Eilthough  Ug  is  the  largest  prime  factor  of  n/n^  .) 

The  average  value  of  a  logarithm  may  seem  at  first  to  be  of  limited 
practical  Interest,  by  comparison  with  the  median  and  other  percentiles; 
however,  we  can  interpret  it  meeuiin  fully  by  saying  that  Dj^m  is  the 
asymptotic  average  number  of  digits  in  the  k-th  largest  prime  factor  of 
an  m-dlglt  number.  Dlckman's  constant  eirises  also  in  an  unexpected 
way  in  connection  with  our  simple  factoring  algorithm:  The  probability 
that  ng  <  Vn^  ,  namely  the  probability  that  the  algorithm  needs  to 
divide  by  all  numbers  up  to  “y/n^  ,  is 


(  f  "iCsiw)  =  {  f  p(^)  =  ( *V(u*i)  (9.8) 


4o 


by  substituting  u  =  2/t  -  1  .  So  this  probability  equals  In 

the  empirical  tests  which  led  to  Figure  2,  exactly  6l  of  the  100 
numbers  had  n^  <  '\/^  . 


Ul 


10.  Relation  to  permutations. 


The  numerical  value  of  in  (9.5)  leads  again  to  a  feeling  of 

deja  vu;  and  sure  enough  Dlckman's  constant  turns  out  to  be  the  same  as 

"Golorab's  constant",  which  has  been  evaluated  to  55  places  in  [6]. 

Golomb's  constant  \  is  defined  to  be  lira  f  /n  ,  where  I  is 

n  -*00  n '  n 

the  average  length  of  the  longest  cycle  in  a  random  permutation.  In 
Golomb's  original  analysis  [5]  of  this  combinatorial  problem  (which  is  not 
obviously  related  to  prime  factors  at  all'.),  he  independently  defined  a 

CO 

_  2 

function  essentially  identical  to  p(a)  ,  and  he  computed  X  =  J  p(t-l)dt/t 

1 

GD 

numerically.  Another  expression  \  =  J  exp(- x  -  E(x) )  dx  was  found  later 

0 

by  L.  Shepp  and  S.  P.  Lloyd  [12]. 

In  Table  1  of  their  paper,  Shepp  and  Lloyd  list  also  the  limiting 

values  i'^'/n  -  f  E(t)^”^  exp(-t -E(t))dt/(k-l) '.  for  the  average  length 
0 

of  the  k-th  longest  cycle;  and  this  agrees  numerically  with  Dj^  for 
1  <  k  <  5  •  In  fact,  the  Shep^  -  Lloyd  formula  yields  Dj^  for  all  k  , 
since 


00 


/ 


(k-l)i 


exp(-t  -E(t))dt 


J  t  e"'‘^(ej^(t)  -  ej^_j^(t))dt 

J  t  e""*^  /  ( Pjj(u)  -  Pij.i(u) )  du  dt 
/  (pjj(u-l)  -pij_i(u-l))  jf  te"'^^'^^  dtdu 

CO 

f  { Pi5.(u-1)  -  Pjj_l(u-1) )  du/u^  .  (10.1) 


U2 


Therefore,  if  we  are  factoring  a  reindom  m-digit  number,  the  distribution 
of  the  n\jmber  of  digits  in  its  prime  factors  is  approximately  the  same  as 
the  distribution  of  the  cycle  lengths  in  a  random  permutation  on  m 
elements'.  (Note  that  there  are  approximately  In  m  factors,  and 
approximately  In  m  cycles.) 

There  is  a  fairly  simple  explanation  for  the  fact  that  Pj^(Q!)  turns 
up  in  the  study  of  cycles  in  permutations.  Let  Qj^(n,r)  be  the  namber 
of  permutations  on  n  objects  having  less  than  k  cycles  of  length 
exceeding  r  .  Then,  by  considering  the  permutations  on  ntl  elements 
(0,1,  ...,n]  and  considering  the  n'./(n-m)'.  possible  cycles  in  which  0 
appears  with  m  different  elements,  we  have 

9  (n.l,r)  =  £  .  (10.2) 

0<m<r  '  '  r<m<n  '  ' 

Therefore  if  qj^(n,r)  =  Qj^(n,r)/n'.  is  the  probability  that  the  k-th 
largest  cycle  has  length  <  r  ,  we  have 

(ttfl)qj^(ntl,r)  =  L  qj^(n-m,r)  +  D  qj^_^(n-ra,r)  ;  (10-5) 

0<m<r  r<m<n 

replacing  n  by  n-1  yields 

nqj^(n,r)  =  E  qj^(n-l-m,r)  +  E  qj^_j_(n-l-m,r)  .  (10.4) 

0<ra<r  r<m<n 

Subtracting  these  two  equations,  we  have 

(n<-l)(qj^(n<-l^r)  -qj^(n,r))  =  qj^_^(n-r,r)  - qj^(n-r,r)  ,  (10. 5) 

and  this  is  analogous  to  the  differential  equation 

45 


a  Pi(0!)  =  Pk.i(«-1)  -  PkCo'-l) 


(10.6) 


The  connection  between  the  two  problems  is  completed  by  showing  that 
=  Pk(n/r) +  0(l/r)  . 

A  similar  distribution  is  obtained  for  the  degrees  of  the  factors 

of  a  random  polynomial  of  degree  n  ,  over  a  finite  field:  The  average 

degree  of  the  k-th  "largest"  irreducible  factor  will  tend  to  be  approximately 

D,  n  . 
k 

Let  us  close  by  stating  an  open  problem:  Are  the  functions 
algebraically  independent?  They  are  linearly  independent,  because  of 
Eq.  (5.5)- 
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Appendix  A.  The  number  of  prime  factors. 


Following  the  notation  of  Hardy  and  Wright  [  ^>  ],  let  (u(n)  be 
the  number  of  distinct  prime  factors  of  n  ,  and  let  n(n)  be  the 
total  number  of  prime  factors  including  multiplicity.  Thus,  n(n) 
is  the  quantity  T  in  the  analysis  of  the  algorithm  above.  Clearly 
1  <  Q(n)  <  logg  n  ,  and  both  of  these  limits  are  obtained  for  infinitely 
many  n  ;  similarly  (u(n)  can  get  as  large  as  In  n  /  In  In  n  .  On 
the  other  hand  these  extreme  values  are  relatively  rare,  and  the  number 
of  factors  is  usually  near  In  In  n  . 

P.  Erdos  and  M.  Kac  [  4  ]  proved  that  the  number  of  n  in  the  range 
1  <  n  <  N  such  that  u)(n)  <  In  In  N+  c  vln  In  N  is 


We  might  say  that  (u(n)  behaves  essentially  like  a  normally  distributed 
random  variable  with  mean  and  variance  In  In  n  ,  where  n  is  large. 

Erdos  and  Kac  remarked  that  their  methods,  which  were  based  on 
the  idea  that  residues  modulo  distinct  primes  are  independent,  could 
be  extended  to  the  case  of  prime  factors  with  multiplicities  included, 
but  they  did  not  state  what  the  resulting  theorem  would  be.  Fortunately 
it  is  easy  to  deduce  the  asymptotic  behavior  of  Q(n)  from  that  of 
u)(n)  ,  using  a  method  like  that  in  [  5  ].  Let  k(N)  be  the  nmber  of  n 
in  1  <  n  <  N  such  that 


Us 


u)(n)  <  1ji  In  N  +  c  Vln  In  N 


(A.5) 


and  let  K(N)  be  the  number  such  that 

n(n)  <  In  In  N  +  c  4ln  In  N  +  In  In  In  N  .  (A. 4) 

Then  lk(N)  -  K(N)  |  is  at  most  the  number  of  n  which  satisiV  (A. 3) 

but  not  (A.1+),  or  (A.U)  but  not  (A.3)>  and  both  of  these  quantities 
are  o(N)  :  If  n  satisfies  (A. 3)  but  not  (A.U),  we  have 
fi(n)  -tt)(n)  >  In  In  In  N  ;  and  the  number  of  such  n  is 
0(n/  In  In  In  N)  ,  because 

E  (Q(n)-u)(n))  =  0(N)  (A.5) 

l<n<N 

by  [6,  Theorem  430].  If  n  satisfies  (A. 4)  but  not  (A. 3)#  then 

In  In  N  +  c  'Jin  In  N  <  m(n)  <  lnlnN  +  [c  +  ^  j  Vln  In  N  , 

\  "Jin  In  N  y 

and  this  is  o(n)  by  the  theorem  of  Erdos  and  Kac. 

We  have  proved  that  the  nvimber  of  n  in  the  range  1  <  n  <  N 
such  that  n(n)  <  In  In  N  +  c  Vln  In  N  is  asymptotically  given  by 
the  normal  distribution  (A.l).  But  this  estimate  is  insensitive  to 
0(1)  terms,  so  the  "average  order"  (6,  Theorem  430]  is  also  relevant: 

lim  ^  E  (u)(n)  -  In  In  N) 

N-®  l<n<N 

=  7  +  E  flogf  1  -  «  -26149  72128  47643  ;  (A.6) 

p  prime  ^  '■ 
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lira  i  Z  (n(n)  -  In  In  N) 

N-»  ^  l<n<N 

=  7+  S  +  l.03l^65  38818  97438  .  (A.?) 

p  prime  ^  ^  P  y  P"  / 

(These  sums  may  be  evaluated  to  high  precision  using  the  formula 

^  S  In  C(ns)  (A. 8) 

p  prime  p  n  >  1 

for  s  >  1  .) 

Let  S  =  [10^*^  -  m  I  1  <  m  <  100}  be  the  numbers  used  to  construct 
Figure  2  above.  For  neS  we  have  In  In  n  «  3*1566  ,  and  the 
following  table  shows  the  actual  distribution  of  u)(n)  and  Q(n)  . 


k  = 

1 

2 

5 

4 

5 

6 

7 

8 

9 

10 

11 

12 

nines  1 

1  «>(n)  =  k)!! 

5 

14 

56 

29 

14 

5 

1 

0 

0 

0 

0 

0 

nines  1 

0 

II 

5 

10 

27 

23 

15 

11 

5 

5 

1 

1 

0 

1 

The  respective  mean  values  are  3*50  and  4.27  •  The  number  of 

square-free  n  (those  with  u)(n)  =  Q(n)  )  was  6l  ,  compared  to  the 

2 

expected  value  600/jt  =  60.795  • 
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Appendix  B.  An  asymptotic  formula  for  . 

In  this  appendix  we  shall  sketch  the  derivation  of  an  asymptotic 
expression  for  Pj(Q!)  as  a  -  <0  .  Our  starting  point  is  the  formula 

a-1  pp(t)dt 

=  [  -Sit- 

.a-1  V  a-l  Po(t)  t^  dt 

'  4  'a'*' I 

0<k<r  of  0  a  0 


we  replace  the  final  term  by  its  asymptotic  value 


,  a-1  (c-t^+ CTt^"'^+ . . .  +  c  ,t)dt  f  -  a-1 

p  1.0 _ 1 _ Ezii_+o  —  r  — 

rH  L  a-t  ^1  Jq  a-t 


0^1  *^0 


>  (B.2) 


so  that  the  remainder  is  0(a  log  a)  .  The  main  integral  in  (B.2) 
is  a  linear  combination  of 


t^  dt  (a-t)^ 

Jq  a-t  =  t 


=  a^(iix  a-Hj^)  -  2  ^  ,  (B.5) 


and  it  remains  to  evaluate  f  pp(t)^  dt  to  0{cr~^  log  a)  . 

0 

Since  pg  =  p^  ,  we  have 


J  P2(t)t^dt 
'0 


J  p(u-l)du  J* 

0  V 


a  k  A 

L 


+  aj^  +  0(a"''"^) 


UQ 


a 


r  p(u-l)(u-l)^  ln(QH-l-u)du 
1 


+  p(u)(u-l)^“'5((Qifl-u)^-l)  +  aj^+0(a‘"‘'^) 

l<?<k  ("’■(j))  Vj/3-* 


l<j<r 


(B.h) 


where  =  J  p(t)  t^dt  =  ACj^  .  Putting  all  this  together  and  summing 


leads  to  the  formula 


2h  2h 

S  (a)  =  (2  lna+l)p  (a) 

a  a 


2b 


+  o(a“^"^)  ,  (b.5) 


a 


where 


In  particular,  <bQ,b^,b2, ...)  =  A^0,2,^,^,^,  > 

"  ...^  .  Since  P5  =  |  (p^Ca)  +  pgC®)  +  Sg(a))  ,  we  have 
the  final  formula 


(B.6) 


p  (a)  =  (lna+l)p  (a) - 

a 


or 


(B.7) 
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